#include  <stdio.h>
#include  <math.h>
#include  <stdlib.h>
#include <sys/io.h>
#include <unistd.h>
#include <string.h>
#include <stdbool.h>
//#include <syslib.h>
#define Npix  2   // 2.8 arcsec
#define Notpix 40  //60, 100*1.7 arcsec. For mini-GWAC, it corresponds to one pixel  error
#define CCDsize_min 100  // 100 for 4k ccd and 2k ccd,  50 for 1k CCD
#define CCDsize_max 4000  // 950 for 1k CCD, 4000 for 4k CCD, 2000 for 2k CCD 
#define datanum 50000
#define datanumnew 50000


main()
{
    int     i,j,k,m,N,N1,N2,N3,F;
	float	rxc1,ryc2,rRA3,rDEC4,rB15,rB26,rR17,rR28,rI9;
	float   xc1,yc2,mag3,merr4,elli5,class6,fwhm7,smag;
    float   otxc1,otyc2,otRA3,otDEC4;
    char    ottype5[200];
    int     flag8;
    float   delx,dely,delr,delotx,deloty,delotr;
    float   delmag;
    float   sum_B2,sum_R2,sum_comp,sum_I;
    float   sum,meandelmag;
    int     Nmatch;
    float   average_comp,average_R2,average_B2,average_I;
    float   x[datanumnew], y[datanumnew],mag[datanumnew],merr[datanumnew],elli[datanumnew],class[datanumnew],fwhm[datanumnew];
    int     flag[datanumnew];
    float   ref1[datanum],ref2[datanum],ref3[datanum],ref4[datanum],ref5[datanum],ref6[datanum],ref7[datanum],ref8[datanum],ref9[datanum];
    float   otx[1],oty[1],otRA[1],otDEC[1];
    char    ottype[100][200];
    FILE    *fp1,*fmatch,*fnomatch,*fcomp,*fdeltamag, *fmeandelmag;
	char    *refall,*newimg,*otcoord;

	refall="gaia.txt";
	newimg="newobj.txt";


	i=j=k=m=N=N1=N2=N3=F=0;
    sum_comp=sum_R2=sum_B2=sum_I=0.0;
//    sum_delmag_R2=sum_delmag_B2=0.0;	
//	fnomatch=fopen("newOT.txt","w+");
	fmatch=fopen("matched.txt","w+");
        fmeandelmag=fopen("meandelmag.txt","w+");
	fp1=fopen(newimg,"r");	
	if(fp1)
	{
	        while((fscanf(fp1,"%f %f %f %f %f\n",&xc1,&yc2,&mag3,&merr4, &elli5))!=EOF)
	        {
                x[i]=xc1;
                y[i]=yc2;
                mag[i]=mag3;
		merr[i]=merr4;
		elli[i]=elli5;
		i++;
		    }
		N1=i;
	}
	fclose(fp1);
    	printf("N1=%d\n",N1);
	
	
    //	 refall is for the usno catalog
	fp1=fopen(refall,"r");
        if(fp1)
        {
          while((fscanf(fp1,"%f %f %f\n",&rxc1,&ryc2,&smag))!=EOF)
          {
          ref1[j]=rxc1;
          ref2[j]=ryc2;
	  ref3[j]=smag;
          j++;
          }
          N2=j;
        }
        fclose(fp1);
        printf("N2=%d\n",N2);
	
    Nmatch=0;
    sum=0;
    meandelmag=0;
    for(i=0;i<N1;i++)
	{
	    k=0;
	    for(j=0;j<N2;j++)
		{
			delx=x[i]-ref1[j];
			dely=y[i]-ref2[j];
			delr=sqrt(pow(delx,2)+pow(dely,2));
        	    	if(delr<=Npix) //this obj is in the usno and in the error circular of miniOT
            		{
                		k++;
				Nmatch++;
				delmag=mag[i]-ref3[j];
				sum=sum+delmag;
                                meandelmag=sum/Nmatch;
				fprintf(fmatch,"%.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f %.3f\n",x[i],y[i],ref1[j],ref2[j],delr,merr[i], elli[i],delmag, mag[i], ref3[j]);	
            			fprintf(fmeandelmag,"%.3f %.3f %d\n",meandelmag,sum,Nmatch);
			}
			break;
		}
//		if(k<1)  //this obj is not in the usno and in the error circular of miniOT
//	    	{            
//			fprintf(fnomatch,"%.3f %.3f %.3f %.3f %.3f\n",x[i],y[i],mag[i],merr[i], elli[i]);	
  //      	}
    	}

        fclose(fmeandelmag);	
//	fclose(fnomatch);
	fclose(fmatch);
}
